A/^-body Simulation of Planetesimal Formation Through 
Gravitational Instability of a Dust Layer 

Shugo Michikoshi^, Shu-ichiro Inutsuka^, Eiichiro Kokubo^, and Izumi Furuya^ 

shugoOtap . scphys .kyoto-u .ac.jp, inutsukaOtap . scphys .kyoto-u. ac . jp , 
kokubo@th.nao.ac.jp, and furuya.izumi@ns-sol.co.jp 

ABSTRACT 

We performed N-body simulations of a dust layer without a gas component 
and examined the formation process of planetesimals. We found that the forma- 
tion process of planetesimals can be divided into three stages: the formation of 
non-axisymmetric wake-like structures, the creation of aggregates, and the colli- 
sional growth of the aggregates. Finally, a few large aggregates and many small 
aggregates are formed. The mass of the largest aggregate is larger than the mass 
predicted by the linear perturbation theory. We examined the dependence of 
system parameters on the planetesimal formation. We found that the mass of 
the largest aggregates increase as the size of the computational domain increases. 
However the ratio of the aggregate mass to the total mass Maggr/Mtotai is almost 
constant 0.8 — 0.9. The mass of the largest aggregate increases with the optical 
depth and the Hill radius of particles. 

Subject headings: methods: N-body simulations, instabilities, gravitation, plan- 
ets and satellites: formation 



1. Introduction 

Many Jupiter-mass planets have been found during recent observations of extrasolar 
planets. The existence of extrasolar planets means that the formation of Jupiter-mass planets 
is a common process. According to the standard model of planet formation, the growth of 
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solid planets occurs in two sequential phases. In the first phase, micron-sized dust grains 
grow into kilometer-sized bodies. These kilometer-sized bodies, called planetesimals, are the 
building blocks of terrestrial planets and cores of giant planets. In the second stage, the 
planetesimals continue to grow through inelastic collisions. 

The initial stage of coagulation of dust grains is supposed to continue up to about 
centimeter size. However, the agglomerative growth of dust grains, from cm to km size, is 
not well understood. The inward orbital drift of meter-sized bodies due to gas drag is very 
rapid. They would fall onto the central star from lAU wit hin a timescale of only 10^ years 
(lAdachi. Hayashi. &: Nakazawalll976l : IWeidenschilling 119771 ). Thus, the direct assemblage of 
dust grains in a laminar flow condition seems difficult. 

One of the scenarios for this stage is the gravitational instability (GI) scenario. If 
particles settle into a layer having sufficiently high density and low velocity dispersion, 
the dust l ayer r apidly collapses owing t o its self-gravity, forming km-sized planetesimals 



(Safronov 



1969 



Goldreich fc WardI Il973l : ICoradini. Magni. fc Federicol Il981 



Sekival Il983 



Yamoto fc Sekiyal |200J) . The timescale of this collapse is only on the order of a Kepler 
period in this scenario, and the rapid migration of meter-sized bodies can be avoided. 

However, there is a critical issue in this GI scenario. As dust grains settle toward the 
midplane, the rotational velocity around the midplane increases because of the reduced effect 
of the gas pressure compared to the centrifugal force and the solar gravity. The rotational 
velocity is a function of vertical distance from the midplane, and the resultant shear may 
induce Kelvin-Helmholtz instability. The slightest amount of turbulence in the nebular 
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layer is difficult if turbulence exists. 



Many authors continue to consider various effects of the GI scenario. iSekiyal (119981 ) 
showed that larg e values of the total dust -to-gas mass ratio may p rovide a density cusp 
at the midplane. lYoudin fc Shul (120021) andlYoudin fc Chiangl (120041 ) argued that this cusp 
can be a trigger of GI. iJohansen et al.l (j2006al ) performed the numerical simulations of the 
Kelvin Helmholtz instability. In the saturated turbulence they foun d a highly overdense 
clump s of dust particles related to the streaming instability found by lYoudin fc Goodman 
(120051 ). iMichikoshi fc Inutsukal (120061 ) performed two fluid analyses of Kelvin-Helmholtz 
instability in the dusty gas layer. They conflrmed the previous result, that the Kelvin- 
Helmholtz turbulence must occur before GI. However, they proposed a possible path toward 
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planetesimal formation through GI, if the dust could grow into a sphere of about 5m . 
Although studies have been conducted on the condition of GI or the linear regime, there has 
been little study of the non-linear regime of GI. The purpose of our study is to examine the 
formation process of planetesimals through using N-body simulations. 

From the linear stability analysis of a two-dimensional self-gravitating disk, the c ondition 
for gr avitational instability of the dust layer is given in terms of Toomre's Q value as (IToomre 
1964h 



Q 



a. 



Sin 



<1, 



11 



where is the radial velocity dispersion, is the Keplerian angular velocity, G is gravita- 
tional constant, and is the surface density of the dust layer. The most unstable wavelength 
Amost of the axisymmetric instability mode for Q = 1 is given by 



A 



most 



02 

"0 



(2) 



If a dust layer becomes unstable, Q < 1, it breaks up into fragments. The mass of a 
planetesimal is estimated as 



M^ 



theor 



A. 



(3) 



At radius a = lAU of the mininium-m ass solar nebula model, the surface density of dust 
is Sd — lOg cm~2. iNakagawa et al.l (119861 ) assumed that the gas is laminar during the settling 



phase and examined the settling path and growth of dust particles due to sweeping process. 
They found that the size of dust particles reaches 20cm. The dynamical optical depth is 
not negligible at lAU ( §2.5p . The direct collision between particles plays an important role 
in the evolution of the dust layer. The dynamical behavior of a self-gravitating collisional 
particle system has been studied in the context of the planetary ring. 

In e arlier studies of planetary rings using N-body simulati ons, gravitational force was ne - 
glect e d (|Brahidll977l) . or only the self-gravity was considered (jWisdom and Tremainelll988h . 
Salol Jl995h extended simulations including mutual gravitational force and performed local 
N-body simulations with various ring parameters. He showed that the wakes are created if 
both self-gravity and inelastic collisions are included. In addition, he found that large surface 
density S and small restitution coefficient e lead to gravitational instability and formation 
of non-axisymmetric wake-like structures and aggregates. Ir i the case of plane tesimal for- 
mation, these aggregates might correspond to planetesimals. iTanga et al.l (l2004l ) performed 
N-body simulations of a gravitationally unstable particle disk in the outer solar system. 
They examined the formation of aggregates through gravitational instability. In the case of 
the particle disk without gas, they confirmed that the wake-like structures are formed. In 
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addition, they performed the simulation including the effect of gas drag. The dissipation 
due to gas drag changed the nonlinear process of Gl. They found the formation of large 
virialized cluster. In their simulation, the optical depth is very small, and the particle disk is 
unstable, even initially. If we consider the formation of planetesimals through gravitational 
instability at — lAU, the dynamic is collision-dominated. In addition, the initial velocity 
dispersion must be large owing to turbulent flow. The particle disk is gravitationally stable 
in the initial condition of our simulations. 



Furuyal (l2004l ) performed a set of simulations in the case of large optical depth. She 
adopted the hard-sphere model as a collision model, and examined parameter dependences 
on the size of computational domain and restitution coefficient. She found that the formation 
process of planetesimal is universal to parameters in the case where initial Toomre's Qinit > 2. 
The formation process is divided into three stages: the formation of non-axisymmetric wake- 
like structures, the creation of aggregates, and the collisional growth of the aggregates. 

In this study, we perform the high-resolution local N-body simulations and study the 
gravitational instability of a dust layer. We extend Furuya's simulations by changing various 
parameters and collision models, and examine the formation process of planetesimals. We 
adopt the hard-sphere and soft-sphere models as a collision model. In the hard-sphere 
model, penetration between particles is prohibited. The velocity changes in an instant when 
a collision occurs. In the soft-sphere model, particles can overlap when they are in contact, 
and the collision has duration. We study the dependences on the size of computational 
domain, the Hill radius of the particle, the optical depth, restitution coefficients, and the 
duration of collision. 

In ^ the simulation method, disk model, and initial conditions are described. Our 
results are presented in §31 We explain the formation process of planetesimals through grav- 
itational instability and examine the time evolution of the typical model and its dependence 
on the disk parameters. In ^ we summarize the results. 



2. Methods of Calculation 

2.1. Equation of Motion 

The effect of gas is important if we consider centimeter-sized dust grains. However, it 
is difficult to consider dust particles and gas consistently. As the first step of our study, we 
neglect the effect of gas for the sake of simplicity. In this paper we examine the nature of 
gravitational instability of a particle disk. This situation corresponds to the disk of large 
particles which are not affected by gas. We will study the effect of gas in the near future. 
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Suppose a central star with mass and dust grains with mass rrii. We introduce 
rotating Cartesian coordinates {x, y, z) called Hill coordinates, the a;-axis pointing radially 
outward, the y-axis pointing in the direction of rotation, and the 2;-axis normal to the 
equatorial plane. The origin of coordinates moves on a circular orbit with the semi- major axis 
ao and the Keplerian angular velocity VLq. We assume rrij -C M, \zj\ -C cq, |%|, \zj\ -C 
aof^o, and \zj\ -C oqVLq, where dot denotes time derivatives. We can write the equation 
of motion in non-dimensional forms independent of mass and the heliocentric distance ao if 
we scale the time by ^, the length by Hill radius ha^ = r^, a nd the mass by h^M^, where 



h = (2mp/3M*)^/^, TUp is the characteristic mass of particles ( Nakazawa fc Idalll988l ). We 
obtain Hill's equation: 

d'^Xi _ dyi \^ _ \ 

dt^ - ^ dt + + 2^ rf/'^^' 

3=1, j^i 
3=1 j^i 

d Zj X ^ m,- 



where r^j is the distance between particles i and j. The first terms on the right-hand side 
are the Coriolis force. The second and third terms denote the tidal force and mutual gravity 
between particles, respectively. 

The equation of motion is invariant under the transformation (x, z, Vx, Vy, Vz) ^ {x + 
nLx, y — {3 /2)nL„Q()t + mL,i, z, v^, v,, — ^QonL ^, f z), where m and n are inte gers, and and 



Ly are the size of the computational domain (jWisdom and Tremainelll988l ). The dynamics 



are independent of the choice of origin. The motion of particles is pursued only in the cell 
with periodic boundary conditions. We use a sliding cell method in which the image comes 
into the cell from the opposite boundary when a particle crosses the boundary of the cell. 



2.2. Gravity Calculation 



In the original Wisdom and Tremaine's method, the mutual gravity was ignored or was 
only treated as the vertical force by enhancement of a vertical frequency. In order to in- 
clude the mutual force exactly, we adopt the subregion method introduced by lDaisaka fc Ida 
(I1999I ). Their method has a cut-off of the gravitational force. ISald ( 1l995l ) studied the influ- 
ence of the cut-off length. A larger cut-off length is needed when the aggregates are formed or 
the wake-like structures appear. The cut-off length in this method is similar to the length of 
a simulation region that is sufficiently large to create the wake-like structures or aggregates. 
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We assume an active region and its copies (see Figured]). Inner and outer cells have 
different angular velocities. These cells slide upward and downward with a shear velocity 
2)VLqLx/2. We divide the simulation region into nine subregions. For each subregion, we make 
virtual regions the same size as the original region. The subregion is centered in the virtual 
region. For particles in the subregion, we calculate gravitational forces from all particles in 
its virtual region that contains copied particles. 

The most expensive part of calculation is the calculation of gravitational forces whose 
cost is Oi N"^). We calculate the gravitational force using the special-purpose GRAPE-6 
computer (IMakino et al.ll2003l ). With GRAPE, calculation time is significantly reduced. 



2.3. Collision Model 



2.3.1. Hard- Sphere Model 



A hard-sphere model has been used in previous simulations of planetary ri ngs (jWisdom and Tremaine 



19881 : ISalol Il99ll : iRichardsonI 1 1994 IPaisaka fc Idal Il999l : IPaisaka et al.ll200lh . In the hard- 
sphere model, particles are considered to be hard and penetration between particles is pro- 
hibited. When a collision event occurs, the particle velocity changes in an instant. In our 
simulation, a collision conserves the relative tangential velocity and reduces the absolute 
value of the relative normal velocity by a factor of e. The relative velocity after a collision 
is described by 

vjj = Vij - (1 + e)(nij ■ Vij)nij, (5) 

where Vy is the relative velocity and rijj is a unit vector along rij. The conservation of 
momentum yields 

mjvj + rra^Vj = mjVi + rra^Vj, (6) 

where v( and vj are velocities of the particle i and j after the collision. Thus the velocities 
after the collision are described by 



m,j 



nii + m 



;i + e)(n 



(7) 



'J 



Vi + 



rrii 



rrii + nij 

If a pair overlaps, r^j < ri + rj, and is approaching ■ Vjj < 0, this pair collides, and change 
their velocities by using Equations ([7]) and ([8]). 



Fig. 1. — Schematic illustration of a computational cell (shaded area surrounded by a thick 
line) with nine subregions (divided by dotted line) for calculation of the mutual gravitational 
force. The region surrounded by a thin line corresponds to the virtual region for a subregion, 
which is created to center the subregion and has the same size as a computational cell. 
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2.3.2. Soft-Sphere Model 



I n the soft-sphere model, particles are permitted to overlap when they are in contact 

( lSaldll995l ). This model allows us to easily calculate the simultaneous collisions of many 
particles. The interaction forces are composed of a restoring harmonic force and a viscous 
damping force, 

F{a) = ka + f3a, (9) 

where a is the penetration depth, k is the stiffness of the string, and (3 is the coefficient of 
viscosity. The force acts in the direction joining particle centers. Assuming that there are 
only two particles, and neglecting the gravitational force, the motions of the pair correspond 
to the damped oscillation whose frequency is given as u = \/u1q — l/(2s)^, where s = Th/f3, 
rh is the reduced mass, and ujq = ^/kjm ( Said 119951 ). The co efficierit of r estitution is 
e = exp (— 7r/(2co's)), and the duration of the impact is = l-njio (jDilleylll993l ). 



2.4. Timestep 

The equatio n of motion is integra ted with a leap-frog scheme. We adopt the variable 



timestep used by iDaisaka &: Idal (Il999l ). In the case of the hard-sphere model, the formula 



of timestep is simply given by 



a; 



At = r/i min -— , (10) 
i ai 



where ?7i is a timestep coefficient, ai and ai are the acceleration of particle i and its time- 
derivative, respectively. In the case of the soft-sphere model, the timestep has to be much 
smaller than the impact duration: 



At = min(?7i min — ^, 772 T!,) (11) 
* I ail 

where 772 is a timestep coefficient. 



2.5. Disk Model 



We assume that all particles have the same mass m^. The system parameters are the 
distance from the Sun ao, the length of region L, the number of particles A^, and the mass 
of particles mp(rp,pp), where pp is the material density of dust particles. The dynam ical 
behavior is characterized by only two non-dimensional parameters ( IDaisaka fc Idalll999l ). the 
dynamical optical depth r and the ratio C = 7'H/2rp. In the case where C > 1, the center 
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of the particle is in the Hill sphere of another particle when two particles come into contact 
with each other. The dynamical optical depth is given by 

r= ^^'^ =01o( ^'^ \( VV (12) 

ippVp ' \lOgcm-2y \2gcm-3y V20cmy 

The ratio ( is given by 

C = 105,528 (^V"7-^V'V-^). (13) 

Using the above two parameters, the other parameters can be estimated. The normalized 
critical wavelength is 

Aniost = 12ittC. (14) 

The normalized size of computational domain L, the number of particle N, the normalized 
radius of particle fp, and Toomre's Q are given by 

L = 127^ArC^ (15) 

N = 5767^AV3C^ (16) 

= ^. (17) 

where A is the ratio L/Amost and dx is the initial radial velocity dispersion scaled by r-nflo. 
In order to determine the size of the computational domain, we have to set A. The initial 
velocity dispersion is also a parameter. We determine the initial velocity dispersion from 
the initial Toomre's Qi^a value using Equation (fT8l) . As shown above, if we set r, C, A, and 
Qinit, we can determine the initial state of the system. 



From Equation ( 1161) . if we set C = 105.528 at lAU, the number of particles becomes too 
la rge to handle w" ith th e cur r ently available computers. We have to use smaller (. According 



to 



Ohtsukil (I1993I ) and I Said (Il995l ). the possibility of accretion is determined by the ratio (. 
Because the capture probabihty decreases abruptly for C < 3/2, the ratio ( should be larger 
than 3/2. Therefore we decided to set C = 2 in our simulation. The effect of ( is discussed 

in mx^ 
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2.6. Initial Condition 

In the Hill coordinates, the Keplerian motion has six parameters: the position of the 



guiding center {xq, Vg), inclination i, eccentricity e, and two phases of motion (INakazawa fc Ida 



19881 ). The eccentricity e and inclination i of particles are assume d to obey the Reyleigh 



distribution with the ratio (e^)^/^/(i^)^/^ = 2 (llda &: Makindll992l ). The eccentricity e is 



related to the radial velocity dispersion (e^)^/^ = \pldxh. The velocity dispersion is deter- 
mined by Equation (fTSj) . The position of guiding center (xg,1/g) and the phases of epi cyclic 
and vertical motions are chosen randomly, avoiding overlapping of particles. 

The initial conditions are summarized in Table [T] Models IH and 1 are the fiducial 
models. In model IH the parameters are r = 0.1, C = 2.0, Qmit = 3.0, and A = 6.0. The 
collision model is the hard-sphere model. Model 1 has the same parameters as model IH 
but its collision model is the soft-sphere model, and the duration of collision is 0.02Tk. 

We can calculate the actual length scales. Using minimum mass solar nebular model the 
most unstable wavelength is 3.3 x lO^km. In the case of the fiducial model, the size of the 
computational domain is 2.0 x lO'^km. If the gaseous disk is turbulent, the vertical structure 
of dust particles is determined by the balance between the sedimentation due to gravity and 
diffusion due to turbulence. The thickness is estimated as H^/Hg = a/{QKT{) where a is a 
non-dimensional facto r of the turbulent viscosity, Tf is the stopping time due to the friction 



( iDubruUe et al.lll995l ). We assume that a is very small as a ~ 10~^. Assuming D = 20cm, 
the equilibrium thickness of a dust layer is 5.6 x lO^km. We take the initial thickness of the 
dust layer in our simulation is about a^/^K = 1-5 x lO^m in the case of the fiducial model. 
The thickness is thinner than the realistic value in order to reduce the calculation time. 

The particle used in this simulation is a "super particle" that represents the group of 
small dust particles at lAU, i.e., the size of particles in this simulation is much larger than 
that of the actual particles and their internal density is much lower than that of actual 
particles. We obtain the radius Vp = 5.9 x lO^cm and the density pp = 1.4 x 10~^gcm~^ from 
( = 2.0, and r = 0.1. One super particle represents 1.5 x 10^^ dust particles. We may think 
that some dust particles in the super particles do not collide when the super particles collide. 
In addition, we may also think that the angles and the relative velocities of the collisions 
of the individual particles in a collision of two super particles may vary much. Thus the 
restitution coefficient e used in this simulation does not represent the physical restitution 
coefficient, but corresponds to the coefficient of the dissipation due to collisions between 
super particles. 
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Table 1. Initial conditions for dust disks 
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cient, the optical depth, the ratio rii/2rp, L/Amost, the duration 
of collision, and the initial Toomre's Q, respectively. 
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2.7. Test Calculation 



We perform the simulations with the same parameters as ISald (119951 ) and lDaisaka &: Ida 



( I1999I ) in order to check the accuracy of our code. We calculate the radial velocity dispersion. 
Each simulation lasts 20Tk and we calculate the average and the standard deviation during 
the last 5Tk. Figure [2] shows the time evolution of the velocity dispersion for the hard-sphere 
and soft-sphere models. The velocity dispersion seems to approach a certain limit value. For 
the parameter ( = 0.82, r = 0.1,0.4,0.6, and e = 0.1,0.4,0.6,0.7, the equilibrium states 
exist. Table [2] shows that the velocity di spers i ons o f our calculation with the s oft-sphere and 
hard-sphere models agree with those of ISalol (119951 ) and iDaisaka fc Idal (119991 ). 



The equilibrium value of the soft-sphere model must approach that of the hard-sphere 
model within the limit of the large stiffness of the spring k —>■ 00. We set the parameters 
as C = 0.82, r = 0.1, and e = 0.5. Figure [2] shows that there is little difference between 
the results of the hard-sphere and soft-sphere models if the duration Tg is sufficiently short. 
The difference is prominent if the aggregates are formed (See §3.3.5p . In the case of these 
parameters, aggregates are not formed, thus we find no remarkable difference. 



Table 2. Equilibrium radial velocity dispersion 



e 


r 


A 


Hard 


Soft 


Salo (1995) 


Daisaka & Ida (1999) 


0.50 


0.40 


6.80 


2.99 ±0.27 


3.05 ±0.30 


3.01 ±0.27 


2.91 ±0.26 


0.50 


0.60 


4.40 


5.06 ± 1.13 


5.30 ± 1.18 


5.84 ± 1.42 


5.09 ± 1.14 


0.10 


0.40 


5.12 


4.10 ±0.83 


4.34 ±0.99 


4.39 ±0.67 


4.50 ±0.59 


0.40 


0.40 


5.12 


3.19 ±0.55 


3.30 ±0.54 


3.48 ± 0.43 


4.06 ±0.59 


0.60 


0.40 


5.12 


2.53 ±0.15 


2.47±0.18 


2.61 ±0.18 


2.58 ±0.17 


0.70 


0.40 


5.12 


3.94 ±0.10 
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3.60 ±0.07 


3.67±0.11 
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Fig. 2. — The example of the time evolution of the velocity dispersion of the hard-sphere 
model {solid curve), the soft-sphere models with Tg = 0.00053 {dashed curve), and = 0.027 
{dotted curve). The parameters are e = 0.5, r = 0.1, and A — 27. 
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3. Results 

3.1. Formation Process of Planetesimals 

We perform the simulation with the various models (see Table [T]). We find that the 
formation process of planetesimals is qualitatively the same in all models, except for the 
model where the initial velocity dispersion is small (Qinit < 1), namely, where a dust layer is 
initially unstable. 

First, we explain the overall formation process of planetesimals through GI in the case 
where Qinit > 1- Figure [3] shows snapshots of model 1 at t/Tx = 0.5,1.0,2.0,3.0,6.0, and 
10.0. In this run, the velocity dispersion of particles is initially large enough for the layer to be 
stable (Qinit = 2.9). The velocity dispersion of particles decreases due to inelastic collisions, 
the dust layer becomes thinner and lessens Toomre's Q value. The non-axisymmetric wake- 
like structure develops gradually when the velocity dispersion becomes suffici ently small 



{Q 2). The wak e-like structure has been examined by many authors (e.g. ISald Il995 



Daisaka &: Idalll999l ). The wakes are created if we include both self-gravity and inelastic 
collisions. The dense parts of wakes fragment into aggregates due to GI. The time scale of 
the fragmentation is the Kepler period. The mass of aggregates is discussed in §3.2.11 These 
aggregates correspond to planetesimals. Once the aggregates are formed, they grow rapidly 
through fast mutual coalescences. In summary, the process of planetesimal formation can 
be divided into three stages: the formation of the non-axisymmetric wake-like structure, 
the formation of aggregates, and the rapid collisional growth of aggregates. The formation 
process for the hard-sphere model (model IH) is the same as that for the soft-sphere model. 

The planetesimal formation in the case with Qinit < 1 is different from the case with 
Qinit > 1 as shown in Figure HI The particle distribution imediately becomes non-uniform, 
and the non-axisymmetric wake structure is not formed. The denser regions of the non- 
uniform structure collapse and directly grow up into aggregates. The subsequent evolution 
is the same as that for Qmit > 1- 

The ffow of the gas component may be turbulent in the initial stage of planetes- 
i mal formation because of Kelvin-Helmholtz instability or magnetorotational instability 



(IBalbus fc Hawleylll99ll ). Dust grains are stirred by turbulent gas, thus their velocity disper- 
sion might initially be large. Therefore, we think the case where Qinit > 2 is more realistic. 
In the following section, we examine the case where Qmit > 2. 
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Fig. 3. — The spatial particle distribution in Model 1 ait — 0.5Tk {top left panel), t — I.OTk 
{top middle panel), t = 2.0Tk {top right panel), t = S.OTk {bottom left panel), t — G.OTk 
{bottom middle panel), and t — IO.OTk {bottom right panel). Solid circles denote the Hill 
sphere of aggregates. 
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Fig. 4. — Same as the Figure [3]but for the case of Qmit = 0.42 (model 4). 
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3.2. Time Evolution 

We investigate the time evolution in our fiducial model (models IH and 1). We describe 
the processes common to all models. We show the time evolutions of mass distribution, 
number of aggregates, mass and radius of the largest aggregate, and the velocity dispersion. 
We make a comparison between the hard-sphere and soft-sphere models. 



3.2.1. Mass Distribution of Planetesimals 



We define aggregates by the method used by iFuruyal (12004 ). At first, we detect an 
aggregate of particles whose mutual distances are smaller than the Hill radius of each particle. 
These particles are strongly attracted by mutual gravitational force. However, some particles 
that have sufficiently large velocity cannot stay in the aggregate. We must impose the binding 
condition and remove grav itationally unbound par ticles. The binding condition is determined 
by the Jacobi energy (e.g. iNakazawa fc Idalll988l ): 
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where v is the relative velocity between the aggregate's center of mass and i-i\i particle, x 
and z are the relative position between the aggregate's center of mass and i-th. particle, n^c 
is the number of particles in the aggregate, and thc is the Hill radius of the aggregate. 

Figure O shows the mass distribution for model 1. The aggregate with mass m > 
O.lmtheor is not formed in < t < 0.3Tk. However, we observe the non-axisymmetric 
wake-like structure at t = I.OTk. The wake-like structure is caused by GI. At t = 2Tk 
the aggregates with mass 0.1 — l.OMtheor are formed. The masses of these aggregates are 
on the same order of fragmentation mass estimated by the linear theory. For t > 3Tk the 
number of aggregates decreases and the larger aggregates are formed by the coalescence of 
the aggregates. Finally, a few large aggregates and somewhat smaller aggregates are formed. 
The largest aggregate reaches about 12Mtheor in this run at t = IO.OTk. The time evolution 
of the mass distribution for the hard-sphere model (model IH) is the same as that for the 
soft-sphere model. 



3.2.2. Number of Aggregates 



We cannot find any remarkable difference between the soft-sphere and hard-sphere mod- 
els in the time evolution of the number of aggregates in Figure [61 For t < 2Tk, the number 
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Fig. 5. — The snapshots of the mass distribution for model 1 at t = 0.5Tk {top left panel), 
t = I.OTk {top middle panel), t = 2.0Tk {top right panel), t = 3.0Tk {bottom left panel), 
t = 6.0Tk {bottom middle panel), and t = IO.OTk {bottom right panel). The vertical axis 
denotes the number of aggregates per unit mass and area. 
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of aggregates increases monotonically up to about 1.2 per A^^g^- According to the linear 
theory, a few aggregates are formed per A^^g^.. This result indicates that these aggregates 
correspond to planetesimals predicted in the linear theory. However, the number of aggre- 
gates decreases rapidly after t ~ 2Tk. This decay is caused by the fast coalescence of the 
aggregates. A few large aggregates absorb many other small aggregates, and thus the total 
number of aggregates decreases. 



The time evolution of the mass of the largest and second largest aggregates is shown 
in Figure [61 In the hard-sphere model (model IH), the largest aggregates grow to about 
Smtheor until t = 4Tk. There is little difference between the largest aggregate and the 
second-aggregate for t < 3Tk. However, the mass of the second- largest aggregate decreases 
abruptly at t ~ 3Tk. This decrease is caused by the coalescence of the largest aggregate and 
the second-largest aggregates. A few aggregates absorbed almost mass in the computational 
domain. 

At t ^ 9.6Tk, there is a deep cusp, that is, the mass of the largest aggregate decreases 
abruptly, and restores immediately. This is an artificial phenomenon due to the cluster 
finding method. We see such a rapid decrease and restoration, when two large aggregates 
collide but do not finally merge into a single aggregate. When a collision between large 
aggregates with large relative velocity occurs, the particles of the transient object have high 
relative velocities, so the binding conditions for the particles in the transient object are not 
satisfied. Thus the transient object is not treated as a bound aggregate defined in §3.2.11 
Finally, the hot transient object breaks up into two large aggregates. 

The Hill radius th and physical radius Vc of the largest aggregate are shown in Figure 
[61 The physical radius is defined as 



where is the position of j-th particle, xq is the position of the aggregate's center of mass. 
The Hill radius th is about 5 times larger than the physical radius Vc- 

In the soft-sphere model (Model 1), the mass of the largest aggregate follows an almost 
similar evolution. However, we do not observe the deep cusp observed in the hard-sphere 
model because there happens to be no impact between the largest and the second-largest 
aggregate in this simulation. Note however that the properties of the largest aggregate is 



3.2.3. The Largest Aggregate 




(20) 
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Fig. 6. — The time evolution for the soft-sphere model (model 1) {left panel) and the hard- 
sphere model (model IH) (right panel). The quantities in this figure are the number of 
aggregates, the mass of the largest and second-largest aggregates, the Hill radius [solid curve) 
and physical radius {dotted curve) of the largest aggregate, and the velocity dispersion of 
field particles from top to bottom, respectively. In the bottom panel we show the velocities 
for Q = 2 {dotted curve) and Q — 1 {dashed curve). 
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sensitive to the random seed in the initial condition. We observed deep cusps in the plot 
that corresponds to a collision between the largest and second-largest aggregates in some of 
the other runs, even with the soft-sphere model. 

3.2.4- Velocity Dispersion 

The time evolution of the radial velocity dispersion is shown in Figure [61 The velocity 
dispersion decreases monotonically until t ~ ITk because of inelastic collisions. When the 
velocity dispersion becomes sufficiently small ((5 — 2), gravitational instability occurs and 
the wake- like structure appears. As a result, the gravitational energy is converted to kinetic 
energy, and the velocity dispersion increases (t > ITk). Once large aggregates are formed, 
other small field particles are scattered by the large aggregates and the velocity dispersion 
of field particles becomes large. Finally, the hot thick disk and the large aggregates on 
the midplane are formed. In reaction, dynamical friction from field particles decreases the 
velocity dispersion of the largest aggregate. 

3.3. Parameter Dependence 

We summarize how the results are affected by changing parameters: the size of the 
computational domain A, the restitution coefficient e, the optical depth r, the ratio of Hill 
radius to the physical diameter (, and the duration of collision Tg. We consider the ratio of 
the aggregate mass to the total mass, the mass of the largest aggregate, the time evolution 
of the physical radius, and the radial velocity dispersion of field particles. 

3.3.1. Size of Computational Domain 

We set r = 0.05 and vary the size of computational domain from A = 3 to 12 by 
keeping the other parameters, (, and Tg at fiducial value (models 12, 13, 14, 15, 16). We 
use shallow optical depth to reduce the number of particles. We check the dependence of 
the final mass of the largest aggregate on the size of the computational domain. As shown 
in Figure U\ the final mass and physical radius of the largest aggregate increase as the size 
of the computational domain increases. The ratios of the aggregate mass to the mass for 
A = 3,6, 8, 10, and 12 are almost constant, and reach about 0.9 in each model. The aggregate 
mass is larger for larger calculation size regions because the total mass increases as the size 
of the computational domain increases. The largest aggregate continues to grow until it 
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absorbs most field particles and small aggregates. 

There is little difference in the radial velocity dispersion of field particles before GI. In 
contrast, the final radial velocity dispersion of field particles increases with the size of the 
computational domain. The scattering of the field particles is more efficient for the larger 
aggregate. Therefore, the velocity dispersion increases as the size of the computational 
domain increases. 

We confirmed that the property of the largest aggregate depends on the size of compu- 
tational domain L. The mass of the largest aggregate increases as the size of computational 
domain L increases. The largest aggregate stops growing owing to the depletion of field 
particles. Thus, the mass of the largest aggregate in our simulation is not realistic. We 
expect that the largest aggregate will possibly stop growing if we perform the simulation on 
a sufficiently large computational domain. We will examine this problem in future work. 



3.3.2. Restitution Coefficient 



The restitution coefficient e is varied from e = 0.1 to 0.9 (models 7, 8, 9, 10, and 11). 
We keep the other parameters, r, ^5 and Tg at the fiducial values. The energy dissipation 
due to the inelastic collisions becomes inefficient as the restitution coefficient e increases. If 
we ignore the energy gain and consider the dissipation due to inelastic collision only, the 
damping timescale of the velocity dis persion is described as Tnarnp — c/ (^^0^(1 — e^)), where c 
is a numerical factor of order of unity (IGoldreich fc Tremaindll978l ) . The damping timescales 
are Tdamp/^K — 1.6,1.7,2.1,3.1,8.4 for e = 0.1,0.3,0.5,0.7, and 0.9, respectively, if we set 
c= 1. 

As shown in Figure [H the gravitational instabilitie s occur only for e = . 1,0.3 , and 



0.5. The critical restituti on coefficient is — 0.7 (e.g. IGoldreich fc Tremaind Il978l : ISalo 



I995I : iDaisaka fc Idalll999l ). When the restitution coefficient e > ecr, the dissipation due to 
collisions is too small to reduce the velocity dispersion. The energy gain from the boundary 
is more efficient than the dissipation. In this case, the velocity dispersion increases mono- 
tonically, and the velocity dispersion cannot reach the critical value. Thus no gravitational 
instability occurs. 

On the other hand, when e < the dissipation is efficient and the velocity dispersion 
can decrease gradually. The decay timescale depends on e and becomes shorter as e decreases. 
The velocity dispersions become minimum when Tmin/TK — 1.4, 1.4,2.3 for e = 0.1,0.3,0.5, 
respectively. This dependence of the minimum timescales Tmin on the restitution coefficient 
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Fig. 7. — The time evolution of the ratio of the aggregate mass to the total mass, the mass of 
the largest aggregate, the physical radius of the largest aggregate, and the velocity dispersion 
of field particles for A = 3 {solid curve), A = 6 {dashed curve), A = 8 {short dashed curve), 
A — 10 {dotted curve) and A — 12 {dash-dotted curve). 
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agrees approximately with that of the damping timescale Tdamp estimated by the analytic 
formula. When the velocity dispersion reaches the critical value, the gravitational instability 
occurs. No remarkable differences are found after the gravitational instability occurs. The 
evolution of the size of the largest aggregate, and the ratio of the aggregate mass to the total 
mass is similar in the models for e = 0.1, 0.3, and 0.5. 

3.3.3. Optical Depth 

We study the effect of optical depth r by comparing results with r = 0.05,0.1, and 
0.125 (models 5, 1, and 6). The collisions occur frequently for larger r. In the case of 
frequent collisions, the dissipation of kinetic energy due to inelastic collisions is efficient. 
The damping timescales of the velocity dispersion of field particles are Tdamp/^K — 3.2, 1.6, 
and 1.3 for r = 0.05,0.1, and 0.125. The damping timescale becomes long as the optical 
depth decreases. 

As shown in Figure [9|, the velocity dispersions of field particles become minimum when 
Train/TK — 3.0, 1.5, and 1.0 for r = 0.05, 0.1, and 0.125, respectively. This dependence of the 
timescales T^in on the optical depth r agrees with that of the damping timescale Tdamp- The 
final velocity dispersion is larger for larger r. This is because the difference in the mass of 
the largest aggregates. The final mass and physical radius of the largest aggregate increases 
as optical depth r increases, as shown in Figure [H This means that frequent collisions lead 
to the efficient growth of the largest aggregate. The ratio of the aggregate mass to the total 
mass is larger for larger r. The ratio reaches about 0.9-1.0 for r = 0.05-1.25. If the collisions 
occur frequently, the field particles grow rapidly. 

3.3.4. Hill Radius 

We set ( = 1.75, 2.0, and 3.0 and fix the other parameters in the fiducial model except 
for the optical depth r = 0.05 (models 17, 18, and 19). We test the effect of the Hill radius. 
The mass and radius of the largest aggregate increase as ( increases in Figure dUl The 
Hill radius corresponds to the distance from which the particle effects the other particles 
through gravitational force. The number of particles that one particle can affect through 
mutual gravitational force increases as ( increases. Therefore, the aggregates may absorb 
field particles more efficiently for larger (. The velocity dispersion of field particles is affected 
by the large aggregates. Therefore, the velocity dispersion of field particles increases as ( 
increases. The ratio of the aggregate mass to the total mass increases more rapidly for larger 
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Fig. 8. — The same as Figure[7]but for e = 0.1 {solid curve), e = 0.3 {dashed curve), e = 0. 
{short dashed curve) , e = 0.7 {dotted dashed curve), e = 0.9 {dash-dotted curve). 
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Fig. 9. — The same as Figure [7] but for r = 0.05 [solid curve), t = 0.1 [dashed curve), 
T = 0.125 [dotted curve) . 
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(. From equation ffTSjl . the critical velocity dispersion for GI is proportional to C^. If we start 
the calculation from the same velocity dispersions, GI occurs earlier for larger (. Therefore, 
if we set larger (, the ratio starts increasing earlier. The ratios reach about 0.7, 0.8, and 0.9 
for ( = 1.75, 2.0, and 3.0. The ratio of the aggregate mass to the total mass is larger if ( is 
larger. 

The realistic Hill radius at lAU corresponds to C ~ 100. As stated in the above 
discussion, a larger aggregate may be formed in the case where ( ^ 100. However, in the 
large ( limit, the mass and radius of the largest aggregate must have upper limits because 
the mass in the computational domain is fixed. We will examine the upper limit of the size 
of the largest aggregate in future work. 

3.3.5. Duration of Collision 

We compare the hard-sphere and soft-sphere models and vary the duration of collisions 
for Ts = 0.01, 0.02, and 0.09 (models IH, 1, 2, and 3) in Figure [TTl The other parameters are 
fixed at the fiducial model. The masses of the largest aggregate reach about 5-lOMthcor- We 
do not observe any clear dependence of the mass of the largest aggregate on the duration of 
collision. This may be caused by the fluctuation. The time evolutions of velocity dispersion 
of the field particles are similar for all models. There is little difference in the velocity 
dispersion for the durations of collisions Ts = 0.01, 0.02, and 0.09. The physical radii for the 
duration of collision Ts = 0.09 is smaller than that for the hard-sphere model. For the soft- 
sphere model, the aggregates can contract and reduce their physical radii. If the aggregates 
reduce their radius, the restoring harmonic force of each particle increases. The balance 
between the restoring harmonic force and the internal gravitational force of the aggregate 
determine its physical radius. The restoring harmonic force decreases as the duration of 
collision increases. Therefore the physical radius for the duration of collision Tg = 0.09 is 
smaller than that for the other models. 



4. Summary 

We performed N-body simulations of the dust layers without a gas component using the 
sliding cell method. The gravitational force that has cutoff length is included by using the 
subregion method. We adopted the hard-sphere and soft-sphere models as collision models. 
We considered the region of the disk around ag — lAU where the dynamics are collision- 
dominated and the Hill radius of particles to be much larger than their physical size. We 
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Fig. 10. — The same as Figure [7] but for ( = 1.8 {solid curve), ( = 2.0 {dashed curve), 
( = 3.0 {dotted curve). 
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Fig. 11. — The same as Figure [7] but for the hard-sphere model {solid curve), and the soft- 
sphere model with Tg = 0.01 [dashed curve), Tg = 0.02 [short dashed curve), and Tg = 0.09 
[dotted curve). 
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examined the formation process of planetesimals through the gravitational instabihty and 
the time evolution of the fiducial model and its dependence on the system parameters. 

If a restitution coefficient e is smaller than the critical restitution coefficient — 0.7, 
though the velocity dispersion is large, the velocity dispersion decreases owing to inelastic 
collisions, and gravitational instabihty occurs when Q '2^ 2. We found that the forma- 
tion process is divided into three stages: the formation of the non-axisymmetric wake-like 
structures, the formation of aggregates, and the collisional growth of the aggregates. The 
non-axisymmetric wake-like structures are formed due to GI. Then the dense parts of wakes 
fragment into aggregates. The aggregates rapidly grow by mutual collisions. Finally, a few 
large aggregates absorb other small aggregates and field particles. 

We found the masses of planetesimals to be larger than these estimated by linear theory. 
The mass range of the aggregates is O.lMtheor ~ 1-OMthcor when the GI occurs. The mass 
distribution changes with time because of the rapid growth of aggregates. The mass of 
the largest aggregates depends on the size of the computational domain and is larger than 
lOMtheor in our calculation. 

We studied the dependence of the results on the system parameters, the optical depth r, 
the ratio of Hill radius to the physical diameter ^, the size of computational domain A, the 
duration of collision Tg, and the restitution coefficient e. We found that the mass and radius 
of the largest aggregate increase as r, C, and A increase. However, the ratio of the aggregate 
mass to the total mass is 0.8-1.0 and almost independent of them. The aggregates absorb 
the most field particles. The total mass in the size of computational domain increases as A 
increases. Therefore the mass of the largest aggregate is larger for larger A. If r is large, the 
collision frequency is large, and accumulation is efficient. The mutual gravitational force is 
stronger for larger C,. Therefore, the growth of the aggregates is more efficient if and r are 
larger. The mass of the largest aggregate is independent of e and Tg. The physical size of 
the largest aggregate is smaller for larger Tg since the restoring force is weaker for larger T^. 

In order to make the simulation of a self-gravitating dust layer feasible, we had to adopt 
some simplification and approximation in the present paper. We used the "super particle" in 
the present simulation. The super particles represent the group of small dust particles. We 
assumed that the super particles obey the simple coUision model analogous to that of the 
small dust particles. The restitution coefficient e used in this simulation is not the physical 
restitution coefficient, but corresponds to the rate of the dissipation due to collisions between 
super particles. Thus, the study of the collision law between super particles and its validity 
remains to be done. 

Unfortunately, the simulation of a self-gravitating dust layer using real-scale dust par- 
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tides is beyond the abihty of today's supercomputers. We used smaller ( = 1.75 — 3.0 than 
the realistic ( ^ 100. We checked that the formation process is common for ( = 1.75 — 3.0, 
and investigated its dependence on (. However, it is necessary to test the model with large 
( to confirm whether our results can be applied to a realistic formation scenario or not. We 
will examine this problem in the near future. 

We found that the masses of planetesimals increase with the size of computational 
domain in our calculation. We did not observe the saturation of this mass increase. If the 
saturation does not occur in the actual formation process, the protoplanets can be formed 
directly through GI. We need to perform the simulation with a much larger size of the 
computational domain to discuss the saturation of the planetesimal growth. 

For the sake of simplicity, we neglected the effect of gas in our simulation. However, the 
effect of gas is not negligible if the dust size is smaller than about 4m and is important for 
the scenario of the formation of planetesimals. Many a uthors have investigated the effect of 
gas: the solid particle s drift radially due to gas drag ( Adachi. Hayashi. fc Nakazawa 19761: 



Weidenschilling 



1976 



1993 



Barge &: Sommerial Il995l : iFromang fc Nelson! l2006l : iJohansen et al 



^^I'^f* ' ^^gul^r momentum dissipation helps gravitational collapse (IWard 

YoudiJ2005aEI). gas turb u lence prevents and helps c o ncentrating solids (IWeidenschilling fc Cuzzil 



20063) . and the 



clumps of solid parti cles are formed due to drag instability flYoudin &: GoodmanI 12005 
Johansen et al.ll2006bl ). In addition, the effect of gas must be important for the nonlinear 
process of the gravitational instability. If the gas is laminar and the turbulence is suffi- 
ciently week, the drag force dissipate the kinetic energy of dust particles. This effect helps 
GI. On the other hand, if the turbulence of gas is strong, the dust particles are stirred up 
and the kinetic energy increases. This effect prevents GI. We expect that the criterion and 
nonlinear process of GI depend on the size of the dust particles and the strength of the tur- 
bulence. We will examine the conditions of gravitational instability and formation process 
of planetesimals in a turbulent gas disk in near future. 

This work is supported by the Grant-in-Aid (No.l5740118, 16077202, 16244120) from 
the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan. 
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